## Bueno de Mesquita and Smith 2009 IO

library(foreign)
library(Amelia)

## Load original dataset (two original datasets merged)
bdms2009 <- read.dta("BDMS2009 IO Rep Data.dta")
head(bdms2009)
dim(bdms2009)

## Drop ID vars
bdms2009$dCNAME <- bdms2009$rCNAME <- bdms2009$countryB <- bdms2009$COLyear <- bdms2009$name <- bdms2009$ccode <- NULL

## Drop derived dummy vars
bdms2009$CAN <- bdms2009$UK <-bdms2009$IRE <-bdms2009$NTH <-bdms2009$BEL <-bdms2009$LUX <-bdms2009$FRN <-bdms2009$SWZ <-bdms2009$SPN <-bdms2009$POR <-bdms2009$GFR <-bdms2009$AUS <-bdms2009$ITA <-bdms2009$GRC <-bdms2009$FIN <- bdms2009$SWE <-bdms2009$NOR <-bdms2009$DEN <-bdms2009$JAP <-bdms2009$ROO <-bdms2009$NZ <-bdms2009$contigkb <- NULL

which( colnames(bdms2009)=="_merge" )
bdms2009 <- bdms2009[, -c(56)]

## How many variables? 110: reduction necessary
dim(bdms2009)

## Which variables are in analysis and have missing data?
sum(is.na(bdms2009 $lnGrossAid))/nrow(bdms2009)*100
sum(is.na(bdms2009 $WA))/nrow(bdms2009)*100
sum(is.na(bdms2009 $RA))/nrow(bdms2009)*100
sum(is.na(bdms2009 $LWB))/nrow(bdms2009)*100
sum(is.na(bdms2009 $LWB2))/nrow(bdms2009)*100
sum(is.na(bdms2009 $LwealthB))/nrow(bdms2009)*100
sum(is.na(bdms2009 $LwealthB2))/nrow(bdms2009)*100
sum(is.na(bdms2009 $lnPOPB))/nrow(bdms2009)*100
sum(is.na(bdms2009 $lnPOPB2))/nrow(bdms2009)*100
sum(is.na(bdms2009 $LkgB))/nrow(bdms2009)*100
sum(is.na(bdms2009 $LkgB2))/nrow(bdms2009)*100
sum(is.na(bdms2009 $LLifeExp))/nrow(bdms2009)*100
sum(is.na(bdms2009 $lnDIS))/nrow(bdms2009)*100 
sum(is.na(bdms2009 $ColdWar))/nrow(bdms2009)*100 ## N
sum(is.na(bdms2009 $colony))/nrow(bdms2009)*100 ## N
sum(is.na(bdms2009 $Ltrade))/nrow(bdms2009)*100
sum(is.na(bdms2009 $Ltau_glob))/nrow(bdms2009)*100
sum(is.na(bdms2009 $Ltau2))/nrow(bdms2009)*100
sum(is.na(bdms2009 $lnMultiAid))/nrow(bdms2009)*100
sum(is.na(bdms2009 $US))/nrow(bdms2009)*100 ## N
sum(is.na(bdms2009 $USRA))/nrow(bdms2009)*100 
sum(is.na(bdms2009 $USLWB))/nrow(bdms2009)*100 
sum(is.na(bdms2009 $USLWB2))/nrow(bdms2009)*100 
sum(is.na(bdms2009 $USLwealthB))/nrow(bdms2009)*100 
sum(is.na(bdms2009 $USLwealthB2))/nrow(bdms2009)*100 
sum(is.na(bdms2009 $USlnPOPB))/nrow(bdms2009)*100 
sum(is.na(bdms2009 $USlnPOPB2))/nrow(bdms2009)*100 
sum(is.na(bdms2009 $USLkgB))/nrow(bdms2009)*100 
sum(is.na(bdms2009 $USLkgB2))/nrow(bdms2009)*100 
sum(is.na(bdms2009 $USLLifeExp))/nrow(bdms2009)*100 
sum(is.na(bdms2009 $USColdWar))/nrow(bdms2009)*100 ## N 
sum(is.na(bdms2009 $USlnDIS))/nrow(bdms2009)*100 
sum(is.na(bdms2009 $UScolony))/nrow(bdms2009)*100 ## N
sum(is.na(bdms2009 $USLtrade))/nrow(bdms2009)*100 
sum(is.na(bdms2009 $USLtau_glob))/nrow(bdms2009)*100 
sum(is.na(bdms2009 $USLtau2))/nrow(bdms2009)*100 
sum(is.na(bdms2009 $USlnMultiAid))/nrow(bdms2009)*100 
sum(is.na(bdms2009 $AnyGrossAid))/nrow(bdms2009)*100 ## N
sum(is.na(bdms2009 $dyad))/nrow(bdms2009)*100 ## N
sum(is.na(bdms2009 $year))/nrow(bdms2009)*100 ## N
sum(is.na(bdms2009 $rCCODE))/nrow(bdms2009)*100 ## N
sum(is.na(bdms2009 $LLifeDATA))/nrow(bdms2009)*100 ## N
sum(is.na(bdms2009 $LRB))/nrow(bdms2009)*100 ## N
sum(is.na(bdms2009 $LRB2))/nrow(bdms2009)*100 ## N

analysis <- as.data.frame(cbind(bdms2009$lnGrossAid, bdms2009$WA, bdms2009$RA, bdms2009$LWB, bdms2009$LWB2, bdms2009$LwealthB, bdms2009$LwealthB2, bdms2009$lnPOPB, bdms2009$lnPOPB2, bdms2009$LkgB, bdms2009$LkgB2, bdms2009$LLifeExp, bdms2009$lnDIS, bdms2009$Ltrade, bdms2009$Ltau_glob, bdms2009$Ltau2, bdms2009$lnMultiAid, bdms2009$USRA, bdms2009$USLWB, bdms2009$USLWB2, bdms2009$USLwealthB, bdms2009$USlnPOPB, bdms2009$USlnPOPB2, bdms2009$USLkgB, bdms2009$USlnDIS, bdms2009$USLtrade, bdms2009$USLtau_glob, bdms2009$USLtau2, bdms2009$USlnMultiAid, bdms2009$USLkgB2,  bdms2009$USLLifeExp, bdms2009$LRB, bdms2009$LRB2, bdms2009$USLwealthB2))

missing <-as.data.frame(cbind(as.integer(complete.cases(bdms2009$lnGrossAid)), as.integer(complete.cases(bdms2009$WA)), as.integer(complete.cases(bdms2009$RA)), as.integer(complete.cases(bdms2009$LWB)), as.integer(complete.cases(bdms2009$LWB2)), as.integer(complete.cases(bdms2009$LRB)), as.integer(complete.cases(bdms2009$LRB2)), as.integer(complete.cases(bdms2009$LwealthB)), as.integer(complete.cases(bdms2009$LwealthB2)), as.integer(complete.cases(bdms2009$lnPOPB)), as.integer(complete.cases(bdms2009$lnPOPB2)), as.integer(complete.cases(bdms2009$LkgB)), as.integer(complete.cases(bdms2009$LkgB2)), as.integer(complete.cases(bdms2009$lnDIS)), as.integer(complete.cases(bdms2009$LLifeExp)), as.integer(complete.cases(bdms2009$Ltrade)), as.integer(complete.cases(bdms2009$Ltau_glob)), as.integer(complete.cases(bdms2009$Ltau2)), as.integer(complete.cases(bdms2009$lnMultiAid)), as.integer(complete.cases(bdms2009$USRA)), as.integer(complete.cases(bdms2009$USLLifeExp)), as.integer(complete.cases(bdms2009$USLkgB2)),  as.integer(complete.cases(bdms2009$USLWB)), as.integer(complete.cases(bdms2009$USLWB2)), as.integer(complete.cases(bdms2009$USLwealthB)), as.integer(complete.cases(bdms2009$USlnPOPB)), as.integer(complete.cases(bdms2009$USlnPOPB2)), as.integer(complete.cases(bdms2009$USLkgB)), as.integer(complete.cases(bdms2009$USlnDIS)), as.integer(complete.cases(bdms2009$USLtrade)), as.integer(complete.cases(bdms2009$USLtau_glob)), as.integer(complete.cases(bdms2009$USLtau2)), as.integer(complete.cases(bdms2009$USlnMultiAid)), as.integer(complete.cases(bdms2009$USLwealthB2))))

dim(analysis)
dim(missing)
apply(missing, 2, sd)

## Remove analysis variables
## bdms2009$lnGrossAid <- bdms2009$WA <-bdms2009$RA <- bdms2009$LWB <- bdms2009$LWB2 <-bdms2009$LwealthB <-bdms2009$LwealthB2 <- bdms2009$lnPOPB <- bdms2009$lnPOPB2 <- bdms2009$LkgB <- bdms2009$LkgB2 <- bdms2009$LLifeExp <- bdms2009$lnDIS <- bdms2009$Ltrade <-bdms2009$Ltau_glob <- bdms2009$Ltau2 <-bdms2009$lnMultiAid <- bdms2009$USRA <- bdms2009$USLWB <- bdms2009$USLWB2 <- bdms2009$USLwealthB <-bdms2009$USlnPOPB <- bdms2009$USlnPOPB2 <- bdms2009$USLkgB <- bdms2009$USlnDIS <- bdms2009$USLtrade <- bdms2009$USLtau_glob <- bdms2009$USLtau2 <-bdms2009$USlnMultiAid <- bdms2009$USLkgB2 <-  bdms2009$USLLifeExp <-  bdms2009$LRB <-  bdms2009$LRB2 <- NULL
## head(bdms2009)
## var(bdms2009)

## Check correlations and missing values
round(cor(bdms2009 $total_g, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $total_g, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $total_g))/nrow(bdms2009)*100
## N
bdms2009$total_g <- NULL

round(cor(bdms2009 $total_n, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $total_n, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $total_n))/nrow(bdms2009)*100
## N
bdms2009$total_n <- NULL

round(cor(bdms2009 $loans_g, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $loans_g, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $loans_g))/nrow(bdms2009)*100
## N
bdms2009$loans_g <- NULL

round(cor(bdms2009 $loans_n, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $loans_n, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $loans_n))/nrow(bdms2009)*100
## N
bdms2009$loans_n <- NULL

round(cor(bdms2009 $grants, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $grants, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $grants))/nrow(bdms2009)*100
## N
bdms2009$grants <- NULL

round(cor(bdms2009 $comm, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $comm, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $comm))/nrow(bdms2009)*100
## N
bdms2009$comm <- NULL

round(cor(bdms2009 $co_gr, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $co_gr, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $co_gr))/nrow(bdms2009)*100
## N
bdms2009$co_gr <- NULL

round(cor(bdms2009 $co_lo, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $co_lo, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $co_lo))/nrow(bdms2009)*100
## N
bdms2009$co_lo <- NULL

round(cor(bdms2009 $food, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $food, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $food))/nrow(bdms2009)*100
## N
bdms2009$food <- NULL

round(cor(bdms2009 $TOTgross, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $TOTgross, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $TOTgross))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $popB, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $popB, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $popB))/nrow(bdms2009)*100
## N
bdms2009$popB <- NULL

round(cor(bdms2009 $kgB, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $kgB, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $kgB))/nrow(bdms2009)*100
## N
bdms2009$kgB <- NULL

round(cor(bdms2009 $rgdpchB, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $rgdpchB, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $rgdpchB))/nrow(bdms2009)*100
## N
bdms2009$rgdpchB <- NULL

round(cor(bdms2009 $GDPB, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $GDPB, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $GDPB))/nrow(bdms2009)*100
## N
bdms2009$GDPB <- NULL

round(cor(bdms2009 $WORLDgdp, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $WORLDgdp, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $WORLDgdp))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $donorGDP, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $donorGDP, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $donorGDP))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $lnGDPB, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $lnGDPB, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $lnGDPB))/nrow(bdms2009)*100
## N
bdms2009$lnGDPB <- NULL

round(cor(bdms2009 $RB, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $RB, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $RB))/nrow(bdms2009)*100
## N
bdms2009$RB <- NULL

round(cor(bdms2009 $regyrB, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $regyrB, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $regyrB))/nrow(bdms2009)*100
## N
bdms2009$regyrB <- NULL

round(cor(bdms2009 $SB, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $SB, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $SB))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $WB, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $WB, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $WB))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $MultiAid, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $MultiAid, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $MultiAid))/nrow(bdms2009)*100
## N
bdms2009$MultiAid <- NULL

round(cor(bdms2009 $popA, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $popA, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $popA))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $kgA, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $kgA, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $kgA))/nrow(bdms2009)*100
## N
bdms2009$kgA <- NULL

round(cor(bdms2009 $rgdpchA, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $rgdpchA, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $rgdpchA))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $GDPA, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $GDPA, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $GDPA))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $lnGDPA, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $lnGDPA, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $lnGDPA))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $lnPOPA, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $lnPOPA, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $lnPOPA))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $A_world, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $A_world, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $A_world))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $A_donor, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $A_donor, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $A_donor))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $regyr, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $regyr, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $regyr))/nrow(bdms2009)*100
## N
bdms2009$regyr <- NULL

round(cor(bdms2009 $SA, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $SA, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $SA))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $tau_glob, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $tau_glob, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $tau_glob))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $s_wt_glo, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $s_wt_glo, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $s_wt_glo))/nrow(bdms2009)*100
## N
bdms2009$s_wt_glo <- NULL

round(cor(bdms2009 $alliance, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $alliance, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $alliance))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $trade, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $trade, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $trade))/nrow(bdms2009)*100
## N
bdms2009$trade <- NULL

round(cor(bdms2009 $logdstab, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $logdstab, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $logdstab))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $McGSmithTrade, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $McGSmithTrade, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $McGSmithTrade))/nrow(bdms2009)*100
## N
bdms2009$McGSmithTrade <- NULL

round(cor(bdms2009 $ipd, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $ipd, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $ipd))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $lntrade, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $lntrade, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $lntrade))/nrow(bdms2009)*100
## N
bdms2009$lntrade <- NULL

round(cor(bdms2009 $Deathrate, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $Deathrate, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $Deathrate))/nrow(bdms2009)*100
## N
bdms2009$Deathrate <- NULL

round(cor(bdms2009 $LDeathrate, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $LDeathrate, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $LDeathrate))/nrow(bdms2009)*100
## N
bdms2009$LDeathrate <- NULL

round(cor(bdms2009 $LSB, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $LSB, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $LSB))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $LRBWB, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $LRBWB, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $LRBWB))/nrow(bdms2009)*100
## N
bdms2009$LRBWB <- NULL

round(cor(bdms2009 $LRB2WB, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $LRB2WB, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $LRB2WB))/nrow(bdms2009)*100
## N
bdms2009$LRB2WB <- NULL

round(cor(bdms2009 $LRBRA, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $LRBRA, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $LRBRA))/nrow(bdms2009)*100
## N
bdms2009$LRBRA <- NULL

round(cor(bdms2009 $LWBRA, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $LWBRA, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $LWBRA))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $LWB2RA, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $LWB2RA, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $LWB2RA))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $LRB2RA, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $LRB2RA, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $LRB2RA))/nrow(bdms2009)*100
## N
bdms2009$LRB2RA <- NULL

round(cor(bdms2009 $USLRBWB, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $USLRBWB, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $USLRBWB))/nrow(bdms2009)*100
## N
bdms2009$USLRBWB <- NULL

round(cor(bdms2009 $USLRB2WB, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $USLRB2WB, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $USLRB2WB))/nrow(bdms2009)*100
## N
bdms2009$USLRB2WB <- NULL

round(cor(bdms2009 $USLRBRA, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $USLRBRA, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $USLRBRA))/nrow(bdms2009)*100
## N
bdms2009$USLRBRA <- NULL

round(cor(bdms2009 $USLWBRA, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $USLWBRA, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $USLWBRA))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $USLWB2RA, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $USLWB2RA, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $USLWB2RA))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $Lbw, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $Lbw, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $Lbw))/nrow(bdms2009)*100
## N
bdms2009$Lbw <- NULL

round(cor(bdms2009 $LRBbw, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $LRBbw, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $LRBbw))/nrow(bdms2009)*100
## N
bdms2009$LRBbw <- NULL

round(cor(bdms2009 $LRB2bw, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $LRB2bw, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $LRB2bw))/nrow(bdms2009)*100
## N
bdms2009$LRB2bw <- NULL

round(cor(bdms2009 $LbwRA, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $LbwRA, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $LbwRA))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $USLbw, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $tot USLbw al_g, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $USLbw))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $USLRBbw, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $USLRBbw, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $USLRBbw))/nrow(bdms2009)*100
## N
bdms2009$USLRBbw <- NULL

round(cor(bdms2009 $USLRB2bw, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $USLRB2bw, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $USLRB2bw))/nrow(bdms2009)*100
## N
bdms2009$USLRB2bw <- NULL

round(cor(bdms2009 $USLbwRA, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $USLbwRA, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $USLbwRA))/nrow(bdms2009)*100
## Y

round(cor(bdms2009 $USLRB, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $USLRB, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $USLRB))/nrow(bdms2009)*100
## N
bdms2009$USLRB <- NULL

round(cor(bdms2009 $USLRB2, analysis, use = "pairwise.complete.obs"), 2)
round(cor(bdms2009 $USLRB2, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(bdms2009 $USLRB2))/nrow(bdms2009)*100
## N
bdms2009$USLRB2 <- NULL

## Imputation
head(bdms2009)
dim(bdms2009)

## What is average percentage of missing data?
NAs <- function(x) {
    as.vector(apply(x, 2, function(x) length(which(is.na(x)))))
    }
NAs(bdms2009)
mean(NAs(bdms2009)/nrow(bdms2009))*100

## Thus: 19 imputations

## Note: lead of LWB already included
set.seed(02138)
bdms2009.out <- amelia(bdms2009, m = 19, idvars = "rCCODE", cs = "dyad", ts = "year", polytime = 3, lags = c("lnGrossAid", "AnyGrossAid", "WA", "RA", "LWB2", "LRB", "LRB2", "lnPOPB", "lnDIS", "LwealthB", "LwealthB2", "LkgB", "LkgB2"), empri = 0.01*nrow(bdms2009))

write.amelia(obj=bdms2009.out, file.stem = "BDMS2009 IO Imp Data", format = "dta", separate = FALSE)